rm(list=ls())
set.seed(8768)
library(fGarch)
library(MASS)
library(nloptr)

#### Output file
nom="Backtest-VaR-exchrate-oos.out"

##### Load the file */
garch = read.csv("exchg.csv",header = FALSE, sep =";")
nmax=ncol(garch)
n=nrow(garch)-1

### log p.d.f.
vrstu<-function(x,nu){ 
if (nu<2){likeli=rep(-1e8,length(x))}
else{lknu=log(gamma((nu+1)/2))-0.5*log(pi*(nu-2))-log(gamma(nu/2));
likeli=lknu-((nu+1)/2)*log(1+(x^2)/(nu-2));
}
return(likeli)}


### Score Standardized student

scorestu<-function(x,nu){
knu=0.5*(-1/(nu-2)+digamma((nu+1)/2)-digamma(nu/2));
fst=(nu+1)/(2*(nu-2))*((x^2)/(nu-2+x^2))-log(1+(x^2)/(nu-2))/2
return(fst+knu)}



    ##Log-likelihood
    garchLLH = function(x,parm) {
        omega = parm[1]; alpha = parm[2]; beta = parm[3];nu=parm[4];
        z =x; ltv = mean(z^2)
        # Use Filter Representation:
        e = omega + alpha * c(ltv, z[-length(x)]^2)
        h = filter(e, beta, "r", init =ltv )
        llh = sum(0.5*log(h)-vrstu(z/sqrt(h),nu))/length(x);
        return(llh)
    }

    #### Log likekihod with constraint  on the parameters
    vrais=function(x,b){
		bb=b;
		bb[1]=b[1]^2;
		bb[3]=exp(-b[2]^2-b[3]^2);
		bb[2]=exp(-b[2]^2)-bb[3];
		return(garchLLH(x,bb))
	}

r=445
maxp=500
stat_exchg=function(num,pp){
u0=garch[,num]
data=log(u0[-1]/u0[1:n])####log-returns
	#Initial conditions
    	data_est=data
	vol=var(data_est)
    	adeb=0.1
    	bdeb=0.8
    	odeb=vol*(1-adeb-bdeb);
	tdeb=c(sqrt(odeb),sqrt(-log(adeb+bdeb)),sqrt(log((adeb+bdeb)/bdeb)),8)

#
######MLE 
    
	log_lik=function(param){vrais(data_est,param)}

    options=list("algorithm"="NLOPT_LN_NELDERMEAD","xtol_rel"=1e-8,"maxeval"=1000)
	res<-nloptr(x0=tdeb,eval_f=log_lik,opts=options)
	estimpar0=res$solution


	


##### Rolling estimator ######
rolling<-function(h){

    data_est=data[h:(r+h-1)]

    #Log likekihod with constraint  on the parameters
    voldeb=var(data_est)
    adeb=0.05
    bdeb=0.9
    odeb=voldeb*(1-adeb-bdeb);
    
    vraisnew=function(x,b,vol){
      bb=b;
      bb[1]=exp(-b[1]^2)*vol
      bb[3]=0.999*exp(-b[3]^2)
      bb[2]=exp(-b[2]^2)*(0.9999-bb[3])+0.00001
      bb[4]=4.01+96*exp(-b[4]^2)
      return(garchLLH(x,bb))
    }
    
    
    k4=mean(data_est^4)/(voldeb^2)
    if (k4>3){nudeb=4.011+6/(k4-3)}else
    {nudeb=4.1}
    if (nudeb>100){nudeb=100}
    tdeb=c(sqrt(-log(1-adeb-bdeb)),sqrt(-log((adeb-0.00001)/(0.9999-bdeb))),sqrt(log(0.999/bdeb)),sqrt(-log((nudeb-4.01)/96)))
    
    ######MLE 
    log_lik2=function(param){vraisnew(data_est,param,voldeb)}
    options=list("algorithm"="NLOPT_LN_NELDERMEAD","xtol_rel"=1e-4,"maxeval"=1000)
    res<-nloptr(x0=tdeb,eval_f=log_lik2,opts=options)
    estimpar=res$solution
    omega=voldeb*exp(-estimpar[1]^2)
    beta=0.999*exp(-estimpar[3]^2)
    alpha=exp(-estimpar[2]^2)*(0.9999-beta)+0.00001
    nuest=4.01+96*exp(-estimpar[4]^2)
    
    nu=nuest

    ###### Calculation of the score function######
	####dlog sigma/dtheta#########
	#### See Zakoian er al. ######
    	z=data[h:(r+h)]
	ltv = mean(z^2)
    	e = omega + alpha * c(ltv, z[-length(z)]^2)
    	hh = filter(e, beta, "r", init =ltv )
    	c3=filter(c(0,hh[-(r+1)]),beta,"r",init=0)
    	c2=filter(c(0,z[-(r+1)]^2),beta,"r",init=0)
    	hh=as.vector(hh)
    	### The derivative of log vol with respect to the parameters.
    	der=cbind(1/hh/(1-beta),c2/hh,c3/hh)
    	der[1,]=c(0,0,0)	

    	######Hit sequence for alpha=p
    	epst=as.vector(as.numeric((z[1:r]/sqrt(hh[1:r]))))
	#### Hits
    	It=1*(pt(epst*sqrt(nu/(nu-2)),nu)<pp)
    	qa=qt(pp,nu)/sqrt(nu/(nu-2))
	coef=qt(pp,nu)*dt(qt(pp,nu),nu)
      ### Correction wrt the true score function
	sc=((nu+1)*(epst^2)/(nu-2+epst^2)-1)
    	fulsc=cbind(sc,scorestu(epst,nu))
    	matV=(t(fulsc)%*%fulsc/r)
    	matV[1,1]=(2*nu)/(nu+3)
    	matP=as.matrix(apply((It-pp)*fulsc,2,mean))
    	matP[1]=-coef #This is MINUS Expectation of "d It/dtheta"
    	dlog=der[1:r,]/2
    	espder=apply(dlog,2,mean)
    	Ahat=cbind(matP[1]*t(espder),matP[2])
    	score=cbind(dlog*sc,scorestu(epst,nu))
    	Vhat=ginv(t(na.omit(score))%*%na.omit(score)/r)
    	corscore=(Ahat)%*%Vhat%*%t(Ahat)
	
	##### Joint moments I_t.I_{t-h} corrected for par. uncertainty
    	Itloc=It-pp	
    	Itm1=c(NA,Itloc[-length(It)])
    	Itm2=c(NA,Itm1[-length(It)])
    	Itm3=c(NA,Itm2[-length(It)])
    	Ahat1=c(mean(Itm1*score[,1],na.rm=TRUE),mean(Itm1*score[,2],na.rm=TRUE),
		  mean(Itm1*score[,3],na.rm=TRUE),mean(Itm1*score[,4],na.rm=TRUE))
    	corscore1=t(Ahat1)%*%Vhat%*%Ahat1

    	Ahat2=c(mean(Itm2*score[,1],na.rm=TRUE),mean(Itm2*score[,2],na.rm=TRUE),
		  mean(Itm2*score[,3],na.rm=TRUE),mean(Itm2*score[,4],na.rm=TRUE))
    	corscore2=t(Ahat2)%*%Vhat%*%Ahat2

    	Ahat3=c(mean(Itm3*score[,1],na.rm=TRUE),mean(Itm3*score[,2],na.rm=TRUE),
		  mean(Itm3*score[,3],na.rm=TRUE),mean(Itm3*score[,4],na.rm=TRUE))
    	corscore3=t(Ahat3)%*%Vhat%*%Ahat3

    #VaR forecasts from R+1
    volf=hh[r+1] # Variance forecast
    epstf=z[r+1]/sqrt(volf) # ratio of real return to vol forecast
    Itf=1*(pt(epstf*sqrt(nu/(nu-2)),nu)<pp)
	############## Our moments ##########
	scf=((nu+1)*(epstf^2)/(nu-2+epstf^2)-1)
    	fulscf=cbind(scf,scorestu(epstf,nu))
	dlogf=der[(r+1),]/2
    	scoref=c(dlogf*scf,scorestu(epstf,nu))

	###### e_t in the paper, i.e. projection onto the orthogonal of an auxiliaray score
    	et=Itf-pp-fulscf%*%ginv(matV)%*%matP
    	# Normalization to have unit variance
    	etf=et/as.numeric(sqrt(pp*(1-pp)-t(matP)%*%ginv(matV)%*%matP))

    	##### etstar in the paper, i.e. projection onto the orthogonal of the true score
    	etstar=Itf-pp-scoref%*%Vhat%*%t(Ahat)
    	# Normalization to have unit variance
    	etstf=etstar/as.numeric(sqrt((pp*(1-pp)-corscore)))

    	##### Correction of the variance of the hit sequence
    	Itbis=(Itf-pp)/as.numeric(sqrt((pp*(1-pp)-corscore)))

    	
	weight<-function(p,r){
		if ((p/r)<1){poids1=-(p/r)^2/3}else{poids1=(2*r)/3/p-1}
		return(poids1)
	}
	calculcov<-function(pp,cor,p,r){
		poids=weight(p,r)
    		return((pp*(1-pp)*pp*(1-pp)+poids*cor)^{0.25})
	}
	calculv<-function(pp,cor,p,r){
		poids=weight(p,r)
    		return((pp*(1-pp)+poids*cor)^{0.5})
	}

	Itcf=Itf-pp
	return(cbind(etf,etstf,Itcf/sqrt(pp*(1-pp)),Itcf/calculv(pp,corscore,500,r),Itcf/calculcov(pp,corscore1,500,r)
,Itcf/calculcov(pp,corscore2,500,r),Itcf/calculcov(pp,corscore3,500,r)))
}
	matrol=t(sapply(seq(1,maxp,1),rolling))

	
	calcultest<-function(vec,nbp){
		x=vec[1:nbp]
		xm1=c(NA,x[-length(x)])
    		xm2=c(NA,xm1[-length(x)])
    		xm3=c(NA,xm2[-length(x)])
      	xix1=(nbp-1)*mean(xm1*x,na.rm=TRUE)^2
    		xix2=(nbp-2)*mean(xm2*x,na.rm=TRUE)^2
    		xix3=(nbp-3)*mean(xm3*x,na.rm=TRUE)^2
		xiM=nbp*mean(x)^2
		return(c(xiM,xix1,xix2,xix3))
	}
	
	##### test P=500
	testet=calcultest(matrol[,1],500)
    	testetstar=calcultest(matrol[,2],500)
    	testItseul=calcultest(matrol[,3],500)
	testIt=500*mean(matrol[1:500,4])^2
	
	x=matrol[1:500,5]
	xm1=c(NA,x[-length(x)])
	testItm1=499*mean(xm1*x,na.rm=TRUE)^2

	x=matrol[1:500,6]
	xm1=c(NA,x[-length(x)])
	xm2=c(NA,xm1[-length(x)])
	testItm2=498*mean(xm2*x,na.rm=TRUE)^2

	x=matrol[1:500,7]
	xm1=c(NA,x[-length(x)])
	xm2=c(NA,xm1[-length(x)])
	xm3=c(NA,xm2[-length(x)])
	testItm3=497*mean(xm3*x,na.rm=TRUE)^2

	test500=c(testet[1],testetstar[1],testIt,testItseul[1],testet[2:4],testetstar[2:4],testItm1,testItm2,testItm3,testItseul[2:4])
	    testos=(1-pchisq(c(test500),1))

return(testos)
    }


outputtable<-matrix(rep(NA,16*(3+nmax*3)),ncol=nmax*3+3)


colnames(outputtable) <-c(" ","UK-US","FF-US","SF-US","Yen-US"," ","UK-US","FF-US","SF-US","Yen-US"," ","UK-US","FF-US","SF-US","Yen-US")
rownames(outputtable) <-c("$e_t$","$e_t^*$","$e_t^c$","I_t","$e_te_{t-1}$","$e_te_{t-2}$","$e_te_{t-3}$",
"$e_t^*e_{t-1}^*$","$e_t^*e_{t-2}^*$","$e_t^*e_{t-3}^*$","$e_t^ce_{t-1}^c$","$e_t^ce_{t-2}^c$","$e_t^ce_{t-3}^c$","$I_tI_{t-1}$","$I_tI_{t-2}$","$I_tI_{t-3}$")


listep=c(0.005,0.01,0.05)
ss=1
is=1
for(pp in listep ){
	ss=ss+1
	is=is+1
	for (numcol in 1:nmax){
	output=stat_exchg(numcol,pp)
	outputtable[,ss]=output[1:16]
	ss=ss+1
}
}
library(xtable)
outputfile=paste(nom,".tex")
print(xtable(outputtable,digits=3),sanitize.text.function=function(x){x},hline.after = getOption("xtable.hline.after", c(4,7,10,13)),file=outputfile)
